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We combine scanning-tunneling-spectroscopy experiments probing magnetic impurities on a su¬ 
perconducting surface with a theoretical analysis of the tunneling processes between (superconduct¬ 
ing) tip and substrate. We show that the current is carried by single-electron tunneling at large 
tip-substrate distances and Andreev reflections at smaller distances. The single-electron current 
requires relaxation processes between the impurity-induced Shiba bound state and the quasiparticle 
continuum, allowing us to extract information on such relaxation processes from our analysis. 


Introduction. —Impurity-induced subgap states pro¬ 
vide a fruitful window into conventional and uncon¬ 
ventional superconductors [1-3]. The Yu-Shiba-Rusinov 
states [4-6] bound by magnetic impurities in conventional 
s-wave superconductors are a simple model system for 
nonmagnetic impurity resonances in unconventional su¬ 
perconductors, probe the competition between supercon¬ 
ducting and Kondo correlations [7-9] , and might provide 
a platform for engineering topological superconducting 
phases with Majorana end states [10-12]. 

In scanning tunneling spectroscopy, Shiba states in¬ 
duce resonances which occur symmetrically at positive 
and negative bias voltages [2, 3, 7, 13]. Given their sub¬ 
gap nature, it is natural to describe the current into Shiba 
states as carried by Andreev processes. These processes 
transfer a Cooper pair into the condensate and are res¬ 
onantly enhanced by the Shiba state [14-17]. Neverthe¬ 
less, STM experiments on Shiba states are typically an¬ 
alyzed in terms of the tunneling density of states which 
is appropriate for single-electron tunneling [1, 18, 19]. 
This allows one to understand the observed asymmetry 
in height between the positive- and negative-bias peaks 
while Andreev processes would necessarily be symmetric 
in bias (for normal-state tips) [17]. 

Here, we combine scanning tunneling mi¬ 
croscopy/spectroscopy (STM/STS) of Shiba states 
using superconducting tips with a comprehensive theo¬ 
retical analysis to elucidate the nature of the tunneling 
processes. We show that both single-electron and 
Andreev tunneling contribute in experiment and explain 
the observed inversion of peak-height asymmetry as 
function of tunneling rates. Our analysis shows that 
STM experiments on Shiba states provide access to 
quasiparticle relaxation rates in superconductors, com¬ 
plementing recent work on superconducting quantum 
dots [20-22] and Josephson junctions [23-26]. 

Experiment .—We have performed STM experiments 
probing Mn adatoms on a Pb(lll) single crystal sur¬ 
face. The experiments were carried out in a Specs JT- 
STM at the base temperature of 1.2 K as well as at 4.8 K. 
The Pb single crystal surface was cleaned by repeated 
sputter/anneal cycles until a clean, atomically flat, and 
superconducting surface was obtained (critical tempera¬ 


ture Tc = 7.2 K and gap A = 1.35 meV at 1.2 K). Mn 
adatoms were evaporated onto the clean sample at a tem¬ 
perature below 10 K, resulting in a density of 30 atoms 
per 100X100 nm^ (see Supplementary Material [27]). Our 
STM experiments were carried out with a Pb-covered, 
superconducting tip (see Ref. [7] for the preparation pro¬ 
cedure) which improves resolution far beyond the Fermi- 
Dirac limit [28, 29]. 

Figure 1 shows spectra of the differential conductance 
dl/dV as a function of bias V, acquired at various 
tip-sample distances and thus tunneling strengths with 
the tip placed above a Mn adatom. All spectra share 
the same characteristic peaks [30] but their intensities 
(normalized to the normal-state conductance) depend 
strongly on the tunneling strength and the sign of the 
bias voltage. 

The peaks in the dl/dV spectra appear at thresh¬ 
olds for various fundamental tunneling processes be¬ 
tween superconducting tip and substrate with mag¬ 
netic adatom: (i) Single electrons can tunnel when the 
negative-energy quasiparticle continuum of the tip over¬ 
laps with the positive-energy continuum of the substrate 
(or vice versa). This requires a threshold voltage eV = 
±2A. (ii) Thermally excited quasiparticles (holes) in the 
positive-(negative-) energy quasiparticle continuum in¬ 
duce a single-particle current even near zero bias, (iii) 
With a Shiba state of energy eq, a single-particle cur¬ 
rent flows when the negative-energy continuum of the 
tip overlaps with the Shiba state, or the positive-energy 
continuum with the symmetric energy — cq. These pro¬ 
cesses have threshold biases eV = ±(A -|- cq)- (iv) Due 
to thermal occupation, a single-electron current can also 
flow when the positive-energy continuum overlaps with 
the Shiba state (and symmetrically when the negative- 
energy continuum overlaps with — cq). This requires a 
threshold bias eV = ±(A — eo). (v) At e\V\ < 2A, an 
electron from, say, the tip can be reflected as a hole, 
transferring a Cooper pair. As all tunneling electrons 
and holes gain an energy eV, (multiple) Andreev pro¬ 
cesses between the quasiparticle continua have thresh¬ 
olds eV = ±2A/n with n = 2,3,.... Andreev processes 
require two or more particles to cross the tunnel bar¬ 
rier and thus become relevant for strong tunnel coupling 
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Figure 1. AI/AV spectra measured on an isolated Mn adatom 
on Pb(lll) for increasing tunneling strength from top to 
bottom (recorded with a lock-in modulation amplitude of 
15 pVjjjjg at a frequency of 912 Hz). Spectra are normalized to 
the ‘normal-state’ conductance measured at 4meV (i.e., well 
outside the superconducting gap) in units of Go = 2e^//i, off¬ 
set for clarity, and scaled when indicated for better visibility. 
The distance to the closest neighboring Mn atom was greater 
than 5nm. A spectrum acquired above the clean Pb(lll) 
surface overlays the smallest-conductance trace (top curve) 
for comparison. The four peaks originating from the deepest 
Shiba level are marked by dashed lines at e\V\ = ±(A ± eo)- 

only [31]. (vi) Shiba states induce additional resonant 
Andreev processes which become relevant at much lower 
tunneling rates. An electron from the negative-energy 
continuum of the tip can virtually tunnel into the Shiba 
state, reflect as a hole, and resonantly transfer a Cooper 
pair into the condensate of the substrate. Together with 
a similar process at reverse bias, this leads to thresholds 
at eV = ±(A -f eo) which coincide with those for the 
single-electron processes. The principal tunneling pro¬ 
cesses involving the Shiba states are sketched in Fig. 2. 

There is an important difference between the single¬ 
electron and resonant Andreev processes [17]. Single¬ 
electron processes change the occupation of the Shiba 
state while Andreev processes merely transfer Cooper 
pairs into the condensate. Thus, a continuous current 
flow by single-electron processes requires relaxation pro¬ 
cesses which empty the Shiba state after it is occupied 
from the tip (or occupy the empty Shiba state), see 
Fig. 2. At finite temperature, a quasiparticle in the Shiba 






Figure 2. Principal tunneling processes involving a Shiba 
state (solid line) within the superconducting gap (enclosed 
by BCS quasiparticle peaks). The chemical potential is rep¬ 
resented by a dashed line, (a) Single-electron tunneling from 
tip to substrate (rate re(a;)) with subsequent relaxation from 
the Shiba state to the quasiparticle continuum (rate Fi). (b) 
Andreev process transfering a Cooper pair to the substrate 
by electron and hole tunneling (with rates re(aj) and rh(uj), 
respectively). The processes in (a) and (b) both contribute 
near the threshold eV = A-|-eo. (c) Single-electron tunneling 
from substrate to tip (with rate re(a;)) after occupation of the 
Shiba state by the relaxation of a thermal quasiparticle (with 
rate r 2 ), contributing to the thermal peak at eV = —(A —eo)- 
The current at the other two thresholds eV = — (A -|- eo) and 
eV = A — eo is carried by analogous hole processes (see [27]). 

state can be excited to the continuum by absorption of a 
phonon or a photon (with rate Fi). Conversely, a ther¬ 
mally excited quasiparticle can relax into the Shiba state 
by emission (with rate r 2 ). 

The observed peaks in the dl/dV spectra can now 
be correlated with Shiba states of energy ~0.22meV, 
~0.77meV, and ~1.18meV, respectively. The multi¬ 
ple Shiba states may reflect different angular-momentum 
channels or spin states S' > 1/2 [13, 32, 33]. To analyze 
the tunneling processes, we focus on the most intense 
Shiba state at Cq — 0.22 meV. This state not only leads 
to the two main peaks at eV = ±(A -|- eg) (with peak 
height a±), but also to two pronounced thermal peaks 
at eV = ±(A — eg) (with peak height j3±). As it is the 
deepest state, its theoretical interpretation turns out to 
be least affected by the presence of the other Shiba states. 

The heights of the peaks associated with this Shiba 
state are plotted in Fig. 3 over several decades in normal- 
state tunneling conductance. We draw attention to two 
important features of these data. First, the peak heights 
vary linearly over a wide region before turning sublin- 
ear at larger tunneling rates. Second, the asymmetry in 
the peak heights a± between positive and negative bi¬ 
ases inverts as a function of tunneling strength: At small 
tunneling rates, a+ < a_, while at large tunneling rates, 
a_|_ > a_. It is also evident that the inversion of the 
peak heights occurs at the crossover between the linear 
and sublinear regimes. 
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Figure 3. Peak heights a± and /3± of the four resonances 
associated with the deepest Shiba level (marked by dashed 
lines in Fig. 1) as a function of normal state conductance 
at T = 1.2K (main panel) and T = 4.8K (inset). The full 
(dashed) lines are fits to Eqs. (1) and (2) for the main (ther¬ 
mal) peaks. The crossover points between single-electron and 
Andreev contributions to a± are indicated by arrows. 


Theoretical analysis .—It is often assumed [1] that the 
peak heights at positive and negative biases measure the 
electron and hole components u and v of the Shiba wave- 
function. The observed inversion of peak heights implies 
that this cannot hold in general. To gain further in¬ 
sight, we calculate the subgap current theoretically by a 
standard Keldysh calculation [34, 35] (see Supplementary 
Material for details [27]). Here, we focus on the physics 
underlying the results. Our calculation includes single¬ 
electron and Andreev processes involving the Shiba state 
as well as phenomenological rates Ti and r 2 for relax¬ 
ation processes between Shiba state and quasiparticle 
continuum. We neglect the non-resonant Andreev re¬ 
flections at the superconducting tip (and thus multiple 
Andreev reflections [36]), which is justified except in the 
regime of very strong tunneling. With this approxima¬ 
tion, the tunneling current becomes a sum of single¬ 
particle and Andreev currents, I = P + 1°“^ with 

_ r duj [ri[renF(w_)-r?,nF(w+)] 

“V 27Th\ (cc - eo)2 + (r/2)2 
r2[re(i - npioj-)) - r/j(i - nF(w+))] 

{OJ - eo )2 + (r/2)2 

duj ThTeinpiuj^) - nF{u!+)] 

2 Trn (w - eo)2 + (r/ 2)2 ' 

Here, the Fermi functions np are evaluated at uj± = w ± 
eV andr = re-i-r;j-i-ri-i-r2. 

The expressions for and P can be understood in 
terms of the basic processes discussed above. The An¬ 




dreev current P involves tunneling of an electron, de¬ 
scribed by re(a;) = 2ttu'^p{lu — eV)t^, and a hole, de¬ 
scribed by r?i(a;) = 2t:v‘^p{oj + eV)t‘^. Here, t is the 
amplitude for tunneling between tip and substrate. The 
rates Fg and F^i are strongly w-dependent through the 
tip’s BCS density of states p{uj). The denominator in 
Eq. (2) reflects the intermediate virtual occupation of the 
Shiba state. It includes the rates Fi for depopulating the 
Shiba state by excitation to the continuum and F 2 for oc¬ 
cupying the Shiba state by a thermally excited quasipar¬ 
ticle. The latter processes are assumed w-independent. 
The four contributions to the single-particle current P 
directly correspond to the peaks [term cx FiFg, see 
Fig. 2(a)], a_ (term cx FiF/j), /?_ [term cx F 2 re, see Fig. 
2(c)], and /3+ (term cx r 2 r?i). 

Eqs. (1) and (2) provide the following basic picture 
consistent with the data in Fig. 3: At weak tunneling, 
the relaxation rates Fi and r 2 are faster than the tip- 
substrate tunneling. Once an electron tunnels into the 
Shiba state from the tip, it is rapidly excited to the 
quasiparticle continuum. In this regime, the tunnel cur¬ 
rent is dominated by the single-electron current 7® which 
is proportional to and thus to the normal-state con¬ 
ductance. The Andreev current P is a small correction 
scaling as This explains the wide linear regime in 
Fig. 3. At stronger tunneling, the tunneling rates be¬ 
come comparable to and eventually larger than the re¬ 
laxation rates F 1 and F 2 . Here, the t-dependence of the 
broadening F leads to a sublinear or even a decreasing 
dependence of the peak heights on the normal-state con¬ 
ductance. As the relaxation processes are thermally ac¬ 
tivated, the crossover point between linear and sublin¬ 
ear regime is strongly temperature dependent, moving to 
lower normal-state conductances for lower temperatures. 
This is consistent with a comparison between main panel 
and inset of Fig. 3. 

Linear regime. —This picture is substantiated by a 
quantitative analysis of the linear regime. For weak tun¬ 
neling and Fi F 2 (i.e., eo ^ T), Eq. (1) yields [27] 

h (ri)3/2 ’ “+ri 

for the peak heights. Here, we introduced the normal- 
state electron (hole) tunneling rate 7 e = 2TTt^VQu'^ {'^h = 
2'Kt^VQV^)^ where vq is the normal-state density of states 
of the tip. The expressions for a_ and /3+ simply differ 
by the substitution u v (or ye -o- y/j). Thus, in this 
regime, the peak height is indeed a measure of the Shiba 
wavefunction at the tip position. From the data in Fig. 3, 
we extract a+/a_ = (rt/n)^ « 0.13. 

All four peaks are related by the relation a+P+ = 
a-P-. This is readily checked against the data in Fig. 3 
and indeed, we find that this identity is well satisfied in 
the linear regime [27]. Moreover, the thermal and main 
peaks in Eq. (3) differ only by a ratio of relaxation rates, 
a+IP- = ri/F 2 = exp(eo/r). Here, the last equality 
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follows from detailed balance. This is in excellent agree¬ 
ment for the data at T = 4.8 K. At T = 1.2 K, we extract 
a slightly higher temperature of T = 1.6 K from the ra¬ 
tio of peak heights. Still, these considerations point to 
a relaxation process involving thermal activation rather 
than the quasiparticle bath suggested in Ref. [17]. 

As r 1 increases with temperature, Eq. (3) also predicts 
the peak heights in the linear regime to decrease with T. 
This is consistent with the data as seen by comparing the 
main panel and the inset in Fig. 3. 

Regime of strong tunneling .—As tunneling rates in¬ 
crease relative to the relaxation rates, the magnitude 
of the single-particle conductance reaches a maximum 
and eventually decreases. As a result, the thermal peaks 
/3± should exhibit a maximum vs normal-state conduc¬ 
tance. The situation is different for the main peaks a± 
with their additional Andreev contribution which keeps 
increasing and eventually dominates the peak magnitude. 
Sufficiently far into this regime, Eqs. (1) and (2) yield 

a+ ~ (2eVM(7/i\/A/eo/(7e’/A)^^^), (4) 

/?_ ~ (2eVM(r2/(7eyA)2/3) (5) 

as well as a_ and /1_|_ which differ again hy u ^ v. The 
main peaks a± keep increasing with tunneling "fe.h, albeit 
with a sublinear dependence. We note that the transition 
from weak to strong tunneling is also accompanied by a 
change in the lineshape of the peaks (see [27]). 

Unlike for normal-metal tips [17, 27], the Andreev 
contribution to the main peaks a± is asymmetric for a 
superconducting tip, but with the asymmetry reversed 
relative to single-electron tunneling. While we have 
a_i_/a_ = {u/vY in the linear regime, Eq. (4) predicts 
a_i_/a_ = in the Andreev-dominated regime. 

Indeed, an inversion of the peak heights a± is seen in 
Fig. 3, as pointed out above. 

Eq. (5) predicts that also the thermal peaks invert, 
from j3-/j5+ = (u/vY in the linear regime to j3-/j5+ = 
in the sublinear regime. This inversion is consis¬ 
tent with the data in Fig. 3. In adddtion, theory predicts 
the thermal peaks to assume a maximum as a function of 
normal state conductance. We observe such a maximum 
only for 13+. For /3_, the peak is expected to occur only 
at rather large normal-state conductance where our ap¬ 
proximations of neglecting multiple Andreev reflections 
and a peak width smaller than eo break down. 

To further substantiate our analysis, we have used Eqs. 
(1) and (2) to fit all four peaks a± and f3± over the en¬ 
tire range of tunneling strengths, see Fig. 3 [27]. There is 
excellent agreement between theory and experiment. We 
attribute the deviations for /3_ at large normal-state con¬ 
ductance to additional contribution from multiple An¬ 
dreev reflections. We can also extract the normal-state 
conductance at which the Andreev and single-particle 
contributions to the main peaks become comparable, see 
the arrows in Fig. 3. (Note that this is distinct from the 


crossover between linear and sublinear dependence.) For 
a+, this happens when 2r;j(2eo) Ti, and for a_, when 
2re(2eo) Ti. As the Andreev contribution 

sets in considerably earlier for a+ than for a_. 

Relaxation rates .—To extract the relaxation rates Ti 
and r 2 quantitatively, we focus on the thermal peak /3+. 
In the sublinear regime it contributes a current I eV 2 /li 
[27]. Moreover, Eq. (3) predicts Ti = {a+/j3-)T2 in the 
linear regime. We can thus extract both relaxation rates 
directly from the experimental data. This yields Ti of 
order 0.2 ns at T = 1.2 K and 6ps at T = 4.8 K. At the 
lower temperature, T 2 differs appreciably from T 1 and we 
extract a value of order 0.7ns. While in principle, one 
could also rely on the main peaks to extract Ti and r 2 , 
this is less accurate since Andreev current and thermal 
peaks also contribute to the total current. 

If the relaxation process relied on directly exciting a 
quasiparticle from the Shiba state to the continuum, we 
would predict a ratio of the relaxation rates at the two 
experimental temperatures of order ~ 10^. We can ac¬ 
count for a substantial part of the apparent discrepancy 
with our observations by recalling that there are addi¬ 
tional Shiba states. If relaxation proceeds as a multistep 
process which first excites to the second Shiba state, we 
predict a ratio of relaxation rates which is consistent with 
experiment (see also [27]). 

Conclusions .—We show that STM experiments on 
subgap states in superconductors probe both single¬ 
electron and Andreev tunneling. We emphasize that such 
experiments are particularly fruitful when performed 
with superconducting tips. In this case, thermal smear¬ 
ing can be neglected and the temperature dependence 
of the current arises entirely from activated quasiparticle 
relaxation processes. Moreover, the additional thermal 
peaks facilitate the analysis and provide access to the re¬ 
laxation rates. We find that at weak tip-substrate tunnel¬ 
ing, the current is dominated by single-electron tunneling 
and linear in the normal-state conductance. This regime 
can be used to map out the bound-state wavefunction. At 
stronger tip-substrate tunneling, the dependence on the 
normal-state conductance becomes sublinear. While the 
dependence on the Shiba wavefunction becomes more in¬ 
volved, this regime provides access to pertinent quasipar¬ 
ticle relaxation rates involving the subgap states. Specifi¬ 
cally, we can extract the rates for quasiparticle relaxation 
into and out of the bound state. The present experiment 
was restricted to two different temperatures. To gain 
further microscopic understanding of the relaxation pro¬ 
cesses, it would be rewarding to perform more systematic 
experiments as a function of temperature. 
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I. EXPERIMENTAL DETAILS 


A. Topography of Mn Adatoms on Pb(lll) 


The Pb(lll) surface was cleaned by Ne’*' sputtering, followed by annealing to T = 430K. This yields an atomically 
clean surface with terraces of several nm width, which are separated by mono-atomic steps [see Fig. S4(a)]. Residual 
Ne atoms from sputtering form nonmagnetic nano-cavities below the surface, which appear as hexagonal protrusions 
or depletions of different sizes in topography [29]. They do not show any signatures of subgap resonances. 

Manganese atoms were deposited on the clean Pb(lll) sample inside the STM at a temperature below 10K. All 
Mn adatoms have the same apparent height after the evaporation. By contacting the adatom with the tip at a bias 
of 5 mV, we can induce a change in the adsorption configuration. The resulting species has a larger apparent height 
[Fig. S4(b) and (c)]. The manipulation is reversible: contact formation at a bias of —180mV results in the initial 
apparent height. After back-manipulation the adatom is shifted laterally with respect to the initial position [see 
line-profiles in Fig. S4(c)]. Thus, the manipulation controllably changes the adsorption site of the adatom. We took 
care that the absolute tip height did not change during the manipulation to rule out a change of the tip apex. 

Both configurations show distinct d//dU spectra. The initial adsorption site has been investigated by Ji et al. [13], 
showing multiple Shiba states. Our spectra on this species show the same characteristic features. In the main 
manuscript we focus on the higher species, because it is stable upon tip approach at low bias and thus allows us to 
investigate dl/dV spectra over a large conductance range. 
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Figure S4. (a) Topography of three terraces of the Pb(lll) surface with an evaporated Mn density of « 30 adatoms per 

100 X 100 nm^. Mn adatoms are marked by black circles. The dark depletions of different sizes are nonmagnetic Neon sub¬ 
surface inclusions, that originate from the sample cleaning process. Setpoint: 50 mV, 200 pA. (b) Topography of a single Mn 
adatom in its two different adsorption states, (c) Line profiles across the atom for all adsorption states. The lateral shift of 
the adatom after the back-manipulation (orange, dotted) is due to jumping into a neighboring adsorption site equivalent to the 
initial one. 


B. d7/dV spectra on a Mn adatom at 4.8 K 

In Fig. S5 we show three examples of dl/dV spectra at the higher temperature of 4.8 K at different tip-sample 
distances, i.e., different tunneling strengths. At this temperature, only two Shiba states are well resolved due to the 
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increased width of the resonances. An additional zero-energy resonance is observed due to tunneling of thermally 
excited quasiparticles in tip and sample [process (ii), as described in the main text]. The peak heights a+, a_, /3+, 
and /?_ as shown in the inset of Fig. 3 were extracted from these spectra. Notice that the relative intensities of a± 
change with increasing tunneling strength. However, unlike at 1.2 K, we do not observe a full inversion of peak heights 
due to the larger relaxation rate. 


W) 

'c 

i— 

3 - 

> 

TD 


-2A -AO A 2A 



Energy (meV) 


Figure S5. AI/AV spectra acquired on a Mn adatom at 4.8K (912Hz, 35p.V,.j„g). Spectra are normalized to the normal-state 
conductance (indicated in the graph). The figure includes assignments of the peaks to the main peaks (a±) and the thermal 
peaks (/3±). 


C. Experiments with a normal-metal tip at 1.2 K 

In the main manuscript, we focus on experiments with a superconducting tip. For completeness, we include spectra 
and corresponding peak-height vs. conductance curves acquired with a normal-metal tip at 1.2 K (Fig. S6). Due to 
thermal broadening of the tip’s Fermi edge (« 360 pV), the energy resolution is drastically decreased compared to 
measurements with superconducting tips and only one pair of Shiba resonances is resolved in Fig. S6(a). From a set 
of such spectra, we extract the peak heights a± [Fig. S6(b)]. In agreement with Ref. [17], we observe an asymmetry in 
the weak coupling regime, which reduces when approaching the Andreev regime and reaches (almost) equal intensity 
for the strongest coupling accessible in the experiment. 


II. THEORETICAL DETAILS 


Here, we derive the expressions for the tunneling current between a superconducting tip and a superconducting 
sample with magnetic impurity, as given in Eqs. (1) and (2) of the main text. We apply the nonequilibrium Green 
function method used in [35]. 
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Energy (meV) normal state dl/dV (Ze^/h) 

Figure S6. (a) Two dl/dV spectra acquired on a Mn adatom at 1.2K with a normal-metal tip (912Hz, 80 Spectra are 

normalized to the normal-state conductance, as indicated in the graph, (b) Peak heights a+ and a_ of the Shiba level as a 
function of normal state conductance at T = 1.2 K, measured with a normal-metal tip. 


A. Green-function expression for the current 


The system is described by the Hamiltonian H = Hl + + Ht, where the three parts describe the tip, the 

substrate, and the tunnel coupling. The superconducting tip is described by the BCS Hamiltonian 


Hl = 


dk 

(27r)3 


,kcr + (a4 



(S6) 


where /2m — p,, /i is the chemical potential, A is the superconducting gap, and Ci^ko- {c/ka) amiiliilates 

(creates) an electron in the tip with momentum k and spin a. The Hamiltonian of the substrate contains a magnetic 
impurity, located at the origin, with spin S pointing along the z direction. The impurity couples to the substrate via 
a potential H5(r) and exchange coupling JSazS{r), where cr^ is a Pauli matrix in spin space. The Hamiltonian takes 
the form 


Hr 


dk 


CfcCfl^ko-Cil.kcr + (Ac^ 

,kt^k-kl 



+ '^{v - JSa)c^j^^^CR^^, 

(7 


(S7) 


where the operator cr^o- = f <ikcR^k<T/(27r)^ annihilates an electron with spin tr at the origin. One can always choose 
a gauge such that the superconducting order parameters in tip and substrate are real. The superconducting phase 
difference then enters the tunneling Hamiltonian 


Ht(t) = [ieA(^)/2c|^^(T)ci^„(T) -b te 

a 


(S8) 


where r is the time argument, t the hopping strength, and we have made the time dependence of Cr/r^„ and „ 
explicit. In writing the tunneling Hamiltonian, we have assumed that the substrate is contacted at the impurity 
location. The time-dependent phase difference between the tip and the sample, 0(r) = 4>o + 2eHr, depends on the 
voltage V applied to the junction. 

The current operator can be obtained from the Heisenberg equation of motion I = —eNj^ = ie[NR,HT\, where Nr 
is the electron-number operator of the tip. We obtain 


/(t) =ie^ \te'-‘^'^"''^/'^c\^{T)cRa(T) - te 

Taking the expectation value yields 

I{t) = eTrjr^ [i{T)G'^^{T, t) - G^j^(r, r)t*(T)] } , 


(S9) 


(SIO) 
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where is a Pauli matrix acting in Nambu space. In the last expression, we introduced the lesser Green function in 
Nambu space 




with i, j = L, R and the hopping matrix 


Ut) = 


^gi 0 (T )/2 Q 

0 


(Sll) 


The time dependence only enters the phase difference and thus the current is a periodic function of time r with period 
2tt/ eV. We can expand the current in a Fourier series in terms of the frequency wq = el^ 




(S12) 


The nonequilibrium Green functions depend on two time arguments and have a generalized Fourier expansion 

G[ri,T2) = / + (S13) 

^ m '' 

We adopt the notation G„m(a;) = G{uj + nu!Q,uj + rnwg) for which the relation G„m(w) = G„_m,o(w + muJo) holds. 
The hopping matrix and its conjugate are given by 


f(r) = 


f*(r) = 


t 0 
0 0 




0 0 

0 -t 


-tUJQT 
^ •) 


0-tJ ^ 


t 0 
0 0 


-tUJQT 
^ •) 


(514) 

(515) 


with inm 


fe/h 

^n—m,0' 


Here, we focus on the dc current which is given by the zeroth order in the Fourier expansion. 


J ^ {^0iGrl,10 01^10 


= h r^' 


,^<,ee 


tG 


<,hh 
RL,-1,0 


-tG 


<,ee 
Li?,01 


-tG 


<.,hh 
LR,0,-1 


(516) 

(517) 


where the superscripts ee and hh denote the two diagonal matrix elements in Nambu space. We do not include the 
nonresonant Andreev reflections at the superconducting tip and thus neglect multiple Andreev reflection processes, 
i.e., = 0, where gL denotes the bare Green function of the tip in the absence of the tunnel coupling. 

Importantly, we retain Andreev reflections at the substrate as they may be resonantly enhanced due to the presence 
of Shiba bound states. We can now write the Green functions Glr appearing in Eq. (S17) in terms of g^ and the 
sample Green functions Gr, which includes tunneling only through the self energy of the Shiba state. Using the 
Langreth rule [34] 


Gl,,L = G^^i*g<+G<i*gl, (SIS) 

G<^ = gliG<+g<iG% (SI9) 


^<,ee ^ ^r,ee 2e ^<ee j_ ^<,ee Te 

'^RL,10 — ^i?,11^10i/L,00 ^i?,11^10yL,00’ 

^<,hh ^ /^r,hh Ih <,hh . ^<.,hh 2h a,hh 

^RL,-1,0 — ^R,-1,-1^-1,o9l,00 ^R,-1,-1^-1,o9l,00^ 

^<,ee ^ r,ee 2e ^<,ee . ^<,ee2e ^a,ee 

^LR,01 — 9l,00^01^R,11 + 9L,00^01^R,in 
^<,hh ^ r,hh2h ^<.,hh . <.,hh2h ^a.hh 

'-'LR,0,-1 — 9l, 00^0,-1'-'R,-l,-l + 9l,oo ^o,-i'-'r,-i,-i- 


(520) 

(521) 

(522) 

(523) 


we obtain 
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Writing G(a; + neV) = Gnn{^), and gL = obtain for the current 

- - G^"^a;_)] g<{u;) - G<’""(a._) [gl{u:) - gl{u;)]} , (S24) 

where we used the short-hand notation uj± = a; ± etd. By using the relation G^ — G^ = G“ — G*", we arrive at 

/ {G>-(a;)g<(a;_) - G<-Hff>(u;_)} 

duj{G>^\u;)g<iu;+) - G<>^\uj)g>{uj+)} . (S25) 


e 


B. Shiba-bound-state Green function 


To determine the Green function Gn of the substrate, we first calculate the bare Green function gji neglecting the 
tunnel coupling to the tip. Without the magnetic impurity, the Green function of a BGS superconductor in Nambu 
space evaluated at the origin is 


5iio(w) = - 


771^0 


\/A 2 - \ ^ ^ J ' 


CO A 


(S26) 


We can include the coupling to the impurity spin in Eq. (S7) by means of the Dyson equation -|- JS — Vtz, 

and obtain 




TTZ/qV A^ — 


(w -I- a-\/A2 — ^2^2 _ ^2 _ _ g_,2^ 

Trr'o 


' -I- (a -|- p)^/ISP- — A 

A w -f (a — j3)\/IsP — uj'^ 


2uja — (1 — -I- 


o; -I- (a -I- j3)y/ISp — uj'^ A 

A w -|- (a — /3)\/A^ — oj"^ ) ' 


(S27) 


where we introduced the dimensionless parameters a = ttuqJS > 0 and /? = Trt'ob^- The subgap states with |a;| < A 
correspond to the poles of the Green function. In particular, in our model the Shiba state energy is given by the pole 

of gn, 


eo = A 


1 - -t /32 


^{1 — 


(S28) 


To calculate the tunneling into the Shiba state, we only need 5 i?,(a;) with uj close to eg. In this limit, we set uj = eg+ Suj 
and expand the denominator in Eq. (S27) to linear order in Scu, 

2uja - (1- a'^ + P'^)ViSP - gp- ~ 2(eo -I- ?HjS)ol - (1 - -I- iSP - Cg ^1 - 


= Slj \ 2a+ 


{l-a^ + P^)eg 


= Suj 


(1 — -I- j3^Y + 4a 


2a 


The numerator can be evaluated at a; = Eq, which leads to 


' + (a±/3)v/A2-aj2:^ A l + (o±/3)" 

V(I-a2+/32)2 + 4c2 


Thus the Green function has the approximate form 


9r{^) = 


1 


a; - eo V MU v 


(S29) 


(S30) 


(S31) 
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with 


2 2 
u ,v = 


2Q;7rt'oA (l + (a ± /l)^) 

((l-a2+/32)2+4^2)3/2' 


(S32) 


Here, u and v are the electron and hole components of the Shiba state (corresponding to the upper and lower sign, 
respectively). Note that in general u ^ v when potential scattering by the impurity is included. 


C. Self energy due to relaxation processes 

Phonon or photon induced relaxation processes introduce a self energy Eph into the substrate Green function 


G = g + gEphff + gEphgEphg + ... 

Approximating the bare substrate Green function g by the contribution of the Shiba state. 


g{u}) = \ips)- 


1 


w - eo 


■(V’sl 


with 


(r|^s) = 


u(r) 
v(r) J ’ 


we find 


G{uj) = IV's)- 


1 


■ii’sl- 


Lu-eo- (V’slSph(w)IV'S') 

We approximate the self energy by its value at a; = eo and retain only the imaginary part, 

Pph = 2Im(V's|Sph(eo)|V’s)- 

Thus, the retarded and advanced Green functions of the Shiba state read 


r,a / \ 

9r = 


1 


a;-eo±irph(eo)/2 \uv v 


u uv 

,2 


(S33) 


(S34) 


(S35) 


(S36) 


(S37) 


(S38) 


Here we have again restricted attention to the Green function at the position of the impurity. 

In quasi-equilibrium, the greater and lesser Green function can be expressed in terms of the retarded and advanced 
Green functions. 


9r{^) = /M(ffK(w) - 9r{^)) = 


S<h(eo) 


(a;-eo)2 + (rph(eo))/2)" V «« 


g>[u:) = -{l-}{u:)){gl-g^^) = 


Sp>h(eo) 


u uv 


(u;-eo)2 + (rph(eo)/2)" V ^ 

where /(w) is the quasi-equilibrium distribution function and we used the relations 

-zs<=rph/, *s> = rph(i -/). 


2 > 


(539) 

(540) 

(541) 


We introduce Pi = zEp^(eo) and r 2 = —zSp^(eo) which can be interpreted as the rates with which the Shiba level is 
emptied or occupied. Note that Pph = Ti -i- r 2 . 
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D. Self energy due to tip-substrate tunneling 


We now include the tunnel coupling of the Shiba state to the tip. For simplicity, we assume that the tip position is 
identical with the impurity position. Then, the self energy due to the tunneling is local at the position of the impurity, 
and we can suppress spatial arguments in the following. The tunneling gives rise to the self energy 


11^10 + ^0-i5£,-1,-1^-10 — 


0 


0 


(S42) 


where we neglect Andreev reflections in the tip as discussed in the main text. Similar relations hold for the self 
energies The retarded and advanced Green functions of the Shiba level coupled to the tip can be obtained 

from the Dyson equation 


f-^r,a _ 1 


-1 r,ax^'r-,a y R , -t^ /o 

^-9r^r u}-eo±iT/2 

where the imaginary part of the self energy leads to a broadening T = re(a;) + r;j(w) + Fi + r 2 with 

re(a;) = 2TTt^u^ p(uj-), 
r,i(a;) = 2T:t^v'^p{uj+) 


1 


,2 

UV V" 


U UV 

,2 


in terms of the BCS density of states 


p{uj) = Vo 


\uj\e{\uj\ - A) 
- A2 


(S43) 


(544) 

(545) 


(S46) 


with vq the normal density of states at the Fermi energy. In Eq. (S43) we have neglected the real part of the self 
energy which would lead to a shift of the resonance energy oc . The lesser Green function of the Shiba state is given 
by [34] 


G<=g<+ + g^^^G^ + 

+ ^rGr) + 9r'^rGr] ■ 


(547) 

(548) 


1 9r'^r 

A straightforward calculation using the relations g^iuj) = 27rip(a;)ni?(a;) and g^iuj) = —2Trip(uj)(l — np^uj)) reveals 

Fa + Fe(a;)np’(w_) + F;j(w)nj?(w+) f 


_ . ^2 1 ^ e\'^\^—) I ^ n\'-^f ti UV 

* (w - eo )2 + (hi + r2 + re(cc) + r^(a;)) 2/4 [uv v\ 


Along the same lines we find 


.Fi + re(w)(l - ni;'(w_)) + Fft,(a;)(l - nF(w+)) f u'^ uv 


* (w-eo)2 + (ri+r2+re(w)+F^H)2/4 Uv v^ 


(S49) 


(S50) 


E. Expressions for the tunneling current 

We can now evaluate the current in Eq. (S25) which yields / = /'’+/“, where 

ef ri[Fe(a;)n_F(a;_) - r,j(a;)nF(w+)] - r2[re(a;)(l - ufU-)) - r,i(a;)(l - nF(w+))] 


nv) = 

nv) = I/d. 


Te{uj)rh{uj) 


(w - eo)2 + (rM/2)' 


(u; - eov + inum 

[nF(w_) - nF(w+)] . 


(551) 

(552) 


The current /“ originates from resonant Andreev reflection, whereas P describes single-particle tunneling and subse¬ 
quent relaxation of quasiparticles in the Shiba state. These two equations are given in the main text as Eqs. (1) and 
( 2 ). 
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F. Relaxation processes 


The intrinsic broadening of the Shiba level is determined by the rates Ti for emptying and r 2 for filling a Shiba 
state. In Ref. [17] the authors considered elastic processes due to a coupling to a fermionic bath as a source of 
quasiparticle relaxation. Another possible relaxation mechanism involves transitions between the Shiba state and the 
quasiparticle continuum assisted by phonons or photons. 

In the present experiment, we find no indication for the presence of a quasiparticle bath at subgap energies as 
the Pb sample exhibits a hard gap away from the impurities. A more likely source are phonon-assisted transitions 
which result in a thermal distribution of the Shiba state / = in the absence of a tunnel coupling, where 

/ = r 2 /(ri -I- r 2 ) [see Eq. (S41)]. Thus for purely thermal relaxation we generally find 


— = 
r2 


(S53) 


The most basic relaxation process involves direct transitions between the Shiba state and the quasiparticle continuum. 
Such processes were studied in Ref. [26] where it was shown that the relaxation rates are given by 


Pi 

r2 



-(A-eo)/T 


(A - eo) -k (A -t eo)e 


-A/T 


(A-eo) + (A + eo)e-^«/^ 


(554) 

(555) 


The relaxation rate Pi for leaving the Shiba state has a thermal factor exp[—(A —eo)/^P] involving the required phonon 
energy of the transition to the continuum A — eg, whereas r 2 is limited by the thermal occupation exp(—A/T) of the 
excited quasiparticles in the continuum. The ratio of the two rates indeed yields Eq. (S53). 

We study the nature of the relaxation processes in the experiment by determining the relaxation rates at different 
temperatures (see discussion in main text and Sec. IV B). While our findings are generally consistent with thermal 
relaxation, the data cannot be readily explained in terms of direct transitions to the continuum. Instead, we propose 
a cascade of thermal relaxation processes via intermediate Shiba states as an alternative mechanism consistent with 
our experimental findings. 


III. CALCULATION OF CURRENT AND DIFFERENTIAL CONDUCTANCE AT THE THRESHOLDS 

In this section, we provide details of the calculations underlying Eqs. (3-5) for the peak conductances in the main 
text. We also calculate the currents at these bias voltages which were used in the main text to extract relaxation 
rates. At the end of this section we illustrate these results by numerically calculating the current and the differential 
conductance as a function of tunneling strength. We organize the calculation by threshold voltages. 


A. eV = A -I- eo 


At this threshold, there are two contributing processes to the current, namely single-electron tunneling into the 
Shiba state as well as resonant Andreev processes: Ia+cq = ^a+eq + -^A-i-eo' 


for the peak current and 


-^A-|-£o 


f(2r^ + ri)^ a;e«ri 

i^(2r;i -|- Pi) We Pl 


a+ 


47re^ 2rh+ri 

h Ti p72 

4-7re^ 2r/,-|-ri 

9h UJe 


We «C Pi 
We 3> Pi 


(S56) 


(S57) 


for the peak conductance. We will now derive these results, first treating the single-electron processes and subsequently 
analyzing the Andreev process. The quantities entering into these expressions will be defined as the calculation 
proceeds. 
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1. Single-electron tunneling 


The relevant single-electron process is shown in Fig. 2(a) in the main text. Analytically, this process contributes 
the term 


lA+eoiV) = el 


duj TiTe{u!)nFi(^ — eV) 
2TTh (w — eo)^-I-(r/2)2 


(S58) 


to the current. Here, the subscript denotes the threshold and the superscript the single-particle (as opposed to 
Andreev) nature of the contributions. The integration variable w denotes the energy of the tunneling electron as 
measured from the Fermi energy of the substrate superconductor. In Eq. (S58) and throughout this section we 
focus on contributions to the current which originate from the vicinity of the BCS singularity. Unlike the differential 
conductance, the full current at the main thresholds e|U| = A -|- eo also includes the contributions from terms oc r 2 
in Eq. (S51), which are responsible for the thermal peaks. The quantitative comparison with experiment in Sec. IV A 
includes all contributions to the current. 

At the threshold eV = A -|- cq, the coupling re(a;) becomes singular exactly at the Shiba energy w = eo because of 
the diverging BCS density of states p{lo — eV) [cf. Eig. 2(a)]. The dominant contribution to the current comes from 
the vicinity of the singularity at cu = eV — A ~ eg and we can approximate 


Teioj) 


2tiv? p{uj 


eV) ~ 7 e 


0(eV-A-uj) 

V 2 y/eV — A — uj 


(S59) 


in terms of the normal-state tunneling rate 7 e = 2TTu^UQt^. In the region of interest the thermal occupation of the 
tip is npiuj — eV) ~ 1. Note that this insensitivity to thermal smearing is a consequence of the superconducting tip. 
(Of course, the current is still sensitive to temperature which enters into the relaxation rates Fi and r 2 .) This also 
implies that the bias voltage enters into the current only via the BCS density of states which is quite distinct from 
the case of a normal-state tip. 

With these ingredients, we can now compute the current in the vicinity of the threshold, 





2Trh 



e{eV-A-uj) 
7 eVA 72 ^ ^ VeU - A - a; ' 

I 


(S60) 


Here, we used that Fi r 2 for eg ^ T and F/j ^ F^. The latter will be justihed below. Note that it is however 
important to keep both Fi and Fg. 

We simplify notation by measuring voltages from the threshold. 


eW = eU-(A + eo), 


and introducing the characteristic energy 



Then, we have 


and 


^A+.o(^') ^ 2ea;3/2p^ f 

Jo 


duj 1 


0 2Trn _ g-p,)2 + (it + 


VZJ ) 


nC 

Jo 


duj 1 


w-eV 


2Trh \/uj 




(S61) 


(S62) 


(S63) 


(S64) 


for the corresponding conductance G = dl/dV. The conductance involves an integral over the tip density of states 
~ multiplied by the function Z{uj) = (uj — eU')/[(a; — eV'Y + (ri/2 + j. We emphasize that Z{uj) 
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Z(w) 


10 
(Oe 

Figure S7. Characteristic function Z{u)) defined in the text at the threshold {eV' = 0) and for Fi = 0. This function determines 
the conductance via Eq. (S64). It is zero at the Shiba state (a; = 0) and peaks at ljs- 


is not the spectral function of the Shiba state as it would be for a normal metal tip. This special feature of the 
superconducting tip arises because the voltage dependence enters through the tip density of states rather than the 
occupation numbers. As shown in Fig. S7 Z(uj) vanishes at the Shiba state (now at w = 0) because of the divergent 
broadening induced by the superconducting tip and also vanishes far from the Shiba energy. In between, it peaks at 
a scale set by the maximum of the effective tunneling strength We and the thermal relaxation rate Fi (in Fig. S7 we 
have set Fi = 0 in which case the peak is at We). 

We can now also discuss the hole tunneling rate r?i(a;) with an associated energy scale ojh = (7h-\/A)^^^/2 in terms 
of the normal state hole tunneling rate jh = 27ru^izot^. In principle, hole tunneling introduces another term into the 
broadening of Z{uj). The broadening then becomes 



3/2 

Fi UJe 
2 


i,ji 


\/w + 2eo 


(S65) 


We can neglect the last term (i.e., the hole contribution to the width of Z{ijj)) as long as max{Fi,a;e} 

In principle, one may imagine situations in which v ^ u so that r/j(a;) contributes significantly to the broadening in 
the strong tunneling regime. As this case is probably irrelevant for this experiment we exclude it from our analytical 
considerations. We discuss implications of a broadening due to r^(tLi) in the presentation of the numerical results at 
the end of this section. 

We focus attention on the peak current and peak conductance. The peak occurs approximately at the threshold 
bias eV = A + eg and we restrict our analytical considerations to the threshold, setting eV' = 0 in the following. While 
this makes our analysis more transparent it also introduces a small numerical error. We emphasize that our results 
exhibit the correct parametric dependence and the quantitative analysis in Sec. IV A is based on the numerically exact 
peak heights. The peak position and height relative to the threshold are discussed in detail in Sec. IIIF. 

Evaluated at the threshold, the integral for /a+eq contains the two energy scales Fi and Wg. For weak tip-substrate 
tunneling, Wg <?C Fi, we can neglect the contribution of Fg to the broadening of Z{uj). In this limit, we find 


^A-I-eo - 


„ 3/2t. duj I 1 

2y/2e Wg/^ dx 1 

ttH Jq + l 


(S66) 


The integral is elementary and we obtain the result 


^A-I-eo - 


2e Wg^^ 

X(r^' 


(S67) 
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The peak conductance in the regime We <C Ti can be calculated along the same lines, 


16\/2e^ w, 


dcj- 


3/2 


fx 


h (ri)3/2 


dx-, 


47re^ 

nr (ri)3/2’ 


(S68) 


where the a;-integration is again elementary. 

In the opposite limit of strong tip-substrate tunneling, Wg 3> Ti, we can neglect the contribution of Ti to the 
broadening of Z^uj). In this limit, we find 


'A-|-£o 


Jo 


duj 1 


1 




0 2TTh y/u + ujUuj 
°° dx 1 


Performing the integral yields 


'xh Jo Vx x'^ + 1/x 


eP, 


“ m ■ 

The peak conductance in the regime We ^ Ti can be calculated along the same lines. 


G 


4e^ 


A-S£o /j 


^Pi / dco 

D 

dx 


OJ 


f°° 


0 -I-wf/w] 

\fx 


Jo 


1/a 


47re^ P1 
9/l We 


(S69) 


(S70) 


(S71) 


2. Andreev contribution 


The Andreev current is given by 


in.,(y) = 2ej 


du! re(w)r/j(w)[nF(w — eV) — np{u! -I- eV)] 

2'Kh (w — eo)^-f (r/2)2 


(S72) 


The Fermi functions can be approximated by np{uj — eV) ~ 1 and np{uj + eV) ~ 0, since w ~ eg and eV ~ A -|- e and 
with We <C eg we can approximate r/i(w) by a constant. In the case eg <C A it simply reads 


r/,(w)~r^ = (S73) 

V ^0 

With these approximations, the integrals become equal to those for the single-particle contribution, with the replace¬ 
ment Pi —> Pft in the numerator and an overall prefactor of two. Note that we can still ignore P/i in the broadening 
of Z{u}) under the assumptions spelled out above. This yields the result summarized in Eqs. (S56) and (S57) above. 


B. eV = -{A + eo) 


At threshold the current is again the sum of single-electron and Andreev processes. The relevant single-electron 
process is shown in Fig. S8. Analytically, this process contributes the term 


IlA-eoiV) = -e I 


duj Fir/,(w)np’(w-I-eF) 
2ph {bj - eg)2 -1- (r/2)2 


(S74) 
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'7' 


Figure S8. Single-particle tunneling processes at (a) eV = —(A-|-eo) and (b) eV = A —eo- At the threshold in (a), an additional 
Andreev process contributes to the current (cf. Fig. 2(b) of the main text). 


to the current. Up to overall signs, this differs from the corresponding process near eV = A -b eo discussed above by 
exchanging the roles of re(w) and r/j(a;). The same exchange characterizes the Andreev contribution. Thus, in effect, 
we can obtain the results for this threshold by interchanging rt o n in the expressions for eV = A -|- cq. This yields 


for the peak current and 


for the peak conductance. 


I-A-eo 


f(2r;, + ri)^ w„<ri 

(2r/j-b Ti) w/j ^ Ti 


a- 


45 ^ 2i\+ri 

h Ti ^372 

47re^ 2rh+ri 

9h ujh 


oJh < Ti 
oJh » Ti 


(S75) 


(S76) 


C. eU = -(A-eo) 


At this thermal threshold, only single-electron processes contribute which are shown in Fig. 2(c) in the main text. 
Analytically, this process is described by 


Noting that 


IU+,,{y) = -e 


du! r 2 re(a;)[l — nF{uj — eV)] 
^ (a; - eo)2 + (r/2)2 


(S77) 


1 — nF{u) — eU) ~ 1 = 1 — nF{e\V\ -b w) ~ 1, 


(S78) 


we see that this differs from the expression for the single-electron current at the threshold eU 
factor r 2 /ri. Thus, we obtain 


for the peak current and 


I-A+to 


fr2^ a;e«ri 
^r2 We » Ti 


/?- 


4irU £2 

h Ti 
47re^ r2 
9h CJe 


We < Ti 
We Ti 


A -b eo merely by a 


(S79) 


(S80) 


for the peak conductance. 
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(a) dIldV{2e^lh) (b) /A+,„(eA//j) 



Figure S9. (a) Differential conductance and (b) current at the threshold eV — A + e vs. normal state conductance. The curves 
are obtained numerically from Eqs. (1) and (2) of the main text. The single-particle (red) and Andreev (orange) contributions 
have maxima at distinct values of normal state dl/dV which separate three regimes with dominating relaxation mechanisms 
Ti, oJe, and (tu^/eo)^^^. While analytical expressions for the asymptotes (dashed lines) in the first two regimes are given by 
Eqs. (S56) and (S57), a similar analysis also yields expressions in regime (iii). The parameters are chosen such that all three 
regimes are visible. We have set Fi = 10“^°, r 2 = 0, v? jvo = 0.001, jvo = 1, and eo = 0.3, where all energies are measured 
in units of A. 


D. eV = A — eo 


This thermal threshold is dominated by the contribution shown in Fig. S8 and given by 

duj T 2 Th{<^)[^ — npioJ + eV)] 


lA-eoiV) = el 


27Tn (w - eo)2 + (r/2)2 


(S81) 


This differs from the thermal threshold at eV = —A + eg by the replacement of re(a;) by r;i(a;). Thus, we obtain the 
current at this threshold by the replacement u v. This yields 


for the peak current and 


Ia-cq 


fr2^ ccth<ri 

^r2 ujh > Ti 


P+ 


9 ^ 3/2 

47re^ T 2 

h Ti r3/2 

4-7re^ r2 
9h uJh 


< Ti 
uJh > Ti 


(S82) 


(S83) 


for the peak conductance. 


E. Discussion and numerical results 

In Fig. S9(a) we plot the differential conductance from the single-particle and Andreev currents at the threshold 
eV = A-l-eo according to Eqs. (S51) and (S52) together with the analytical expression in Eq. (S57). We identify three 
regimes as a function of the normal state conductance Gn ~ (2e^/h)t^r'g which exhibit characteristic power-laws as a 
function of tunneling strength. These regimes can be associated with different dominant broadening mechanisms (from 
weak to strong tunneling): (i) intrinsic relaxation Fi, (ii) electron tunneling We, and (iii) hole tunneling (w^/eg)^/^. 

The crossover between regimes (i) and (ii) occurs at a normal state conductance G^, which can be evaluated from 
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the condition Fi ~ Wg. Equivalently the second crossover point G\ is obtained from We 


Gh- 

2e^ 

voA / 


3/2 

h 


aJ 

5 

G%- 

2e2 

<1 

(eo 

s 3/2 

h 

y6 

VA 

) 


and we find 

(584) 

(585) 


Figure S9(a) shows that the full differential conductance peak a+ = dP/dV + dP/dV consists of the sum of 
two terms that peak at different tunneling strengths and thus typically exhibits two peaks as a function of normal 
state conductance. At the first crossover point the single-particle contribution dP/dV reaches a maximum of order 
2e^//i. This is readily understood from the single-particle current in Eq. (S51), which can be viewed as a resonant 
tunneling process through the Shiba state with rates re(a;) and Fi as depicted in Fig. 2(a) of the main text. When 
the effective electron tunneling rate Wg is equal to Fi the conductance reaches a universal value of the order of the 
conductance quantum. For stronger couplings the conductance decreases and the single-particle current shown in 
Fig. S9(b) saturates to a value determined by the relaxation rate. Indeed, Eq. (S56) yields a current P = eTi/SH 
independent of tunneling strength in this regime. At even stronger coupling the Andreev current P exceeds the 
single-particle contribution and thus the total current, P + P, exhibits a shoulder as a function of tunnel coupling. 
We remind the reader that the current at the main thresholds would have additional contributions from the terms 
oc F 2 Fe/^(a;) in Eq. (S51), which we have excluded from our analytical considerations. We evaluate the full current 
for the quantitative comparison between theory of experiment in Sec. IV A. 

We can estimate the normal state conductance at which the Andreev current becomes the dominant contribution 
to the current at the threshold eV = A -|- eo from Fi ^ F/i. We obtain 


26^ i/qA /eoFf 

h V A3 ■ 


(S86) 


In Fig. 3 of the main text we indicate for the positive and negative main peaks by arrows using the parameters 
given in Sec. IV A. 

At the crossover between regimes (ii) and (iii) the effective electron and hole tunneling rates are equal and the 
Andreev contribution to the conductance becomes resonant and reaches a maximum of order 2e^//i. In regime (iii) 
both contributions to the differential conductance decrease with tunneling strength. This peculiar feature arises 
because of the strong energy dependence of the density of states in the superconducting tip and this regime has 
a sizable extension only when and differ by several orders of magnitude. In the experiment this regime is 
presumably limited to very strong tunneling, Gn > 0.1(2e^//i), where our approach ceases to be valid as the peak 
width becomes of the order of eo- Furthermore the extension of this regime is too narrow to observe a decreasing peak 
height. 

The thermal peak /?_ originates entirely from single-particle tunneling and its peak height simply follows the single¬ 
particle contribution to up to a prefactor F 2 /F 1 in regimes (i) and (ii). The remaining peaks a_ and /3_|_ have 
the same qualitative behavior as a+ and /3_ although with different regime boundaries, which are obtained from 
Eqs. (S84) and (S85) by interchanging u ^ v. The peak heights of all four peaks are shown in Fig. Sll(a) as a 
function of normal state conductance. 


F. Lineshape of the Shiba resonance 

In this section we analyze how the lineshape of the Shiba conductance peak is affected by the tunneling strength 
with important implications for fitting experimental dl/dV traces. In the linear regime the single-particle current in 
Eq. (S58) can be approximated near eV = A -f eo by a convolution of the BCS density of states and a Lorentzian of 
width Fi 


Ain(l^) = ^Treu^t^ 


duj 

2Trh 


p{uj 


eV) 


Ti 

{uj-eof+Tl/A- 


(S87) 


This expression has been used previously [7, 13] to fit experimentally measured Shiba resonances. The intrinsic 
Lorentzian lineshape can be obtained by numerical deconvolution of the data with the BCS density of states of the 
tip. 
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sample bias (mV) 


Figure SIO. (a) Lineshape of the differential conductance given by Eq. (S64) as a function of bias voltage near the main 
Shiba resonance at positive bias. The lineshapes differ qualitatively between the linear (Fi ^ We, blue) and sublinear regimes 
(Fi cje, red). Note that single-particle and Andreev currents give rise to identical lineshapes within the approximations 
of our analysis, (b) Experimental dl/dV traces of the Shiba peak near eV = —(A -|- eo) — —1.57meV in the linear (black) 
and sublinear regime (blue). The spectra are normalized to the peak maximum. In the single-electron regime a clear negative 
differential resistance is visible around —1.65 mV. 


In the sublinear regime where the broadening of the Shiba resonance is determined by the tunnel coupling Fg to 
the tip the current reads 


/subiin(F) 2Treuh^ J ^ ■ (S88) 

The tip density of states now also enters the width of the resonance and this expression does not have the form of a 
convolution. Consequently the lineshape of the Shiba resonance changes qualitatively from linear to sublinear regime 
as shown in Fig. SlO(a), where we plot the voltage dependence of the differential conductance given by Eq. (S64). 
Strikingly, the maximum can occur above or below the threshold eV = A -|- eo depending on the tunneling strength. 
This shift must be accounted for when determining the Shiba state energy by htting experimental lineshapes. In 
addition, a characteristic negative differential conductance dip occurs for Ti We but is absent in the opposite 
regime. The disappearance of this dip in the measured lineshape provides a further indication of the crossover 
between weak and strong tunneling regimes. In the experiment the negative differential conductance dip indeed 
vanishes as tunneling strength is increased from linear to sublinear regime as shown in Fig. SlO(b). The peak position 
does not shift which we attribute to a slight suppression of the tip gap with increasing current. 

Our analytical results in the previous subsections refer to the differential conductance exactly at the threshold 
voltages e|E| = A ± eg. The shift of the maximum away from the threshold yields a somewhat larger peak height. 
Given that the peak heights vary by orders of magnitude in the experiment this deviation of at most 35% only affects 
details but is inessential to our central results. Note that we calculate the actual peak height and not the threshold 
values in the quantitative comparison to the experimental data in Sec. IV A. 


IV. THEORETICAL ANALYSIS OF THE EXPERIMENTAL DATA 


Our analysis implies that the subgap transport provides insight into the population dynamics of the Shiba state, as 
governed by the competition between tunneling and quasiparticle relaxation. We now fit the data against the results 
of our model and extract the quasiparticle lifetime in the Shiba state due to thermal relaxation processes. Besides 
demonstrating the validity of our description, this also yields valuable information about the dominant transport 
mechanisms in experiment as a function of the tunnel coupling between tip and sample. 
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normal state dl/dV (2e%) 


Figure Sll. (a) Differential conductance peaks at the thresholds vs. normal state conductance according to Eqs. (S51) and 
(S52). We have used the parameters mentioned in Sec. IV A. (b) Plot of {a+j3+)/{a-l3-) vs. conductance as extracted from 
Fig. 3 in the main manuscript. The ratio stays constant up to 10“^ Go, where the sub-linear regime sets in. 


A. System parameters and theory fit 


Here we provide details of the theoretical hts to the conductance and current at the Shiba peaks as a function of 
normal state conductance, as shown in Fig. 3 of the main text and Fig. S12(a). Several physical parameters can be 
extracted directly from the measured data without fitting. The Shiba energy can be determined from the location of 
the Shiba peaks as a function of bias voltage. For instance, the separation between the two positive bias peaks a_|_ 
and /3+ is 2eo. The same holds for a- and /?_. From the data we estimate cq — 0.22meV for the lowest Shiba level. 
According to Eqs. {S57) and (S76) valid in the linear regime we furthermore obtain (u/v)^ = a+/a_ ~ 0.13. Using 
Eq. (S83) in addition yields ri/r 2 = Q!+//3_ ~ 4. Finally, we can accurately determine r 2 = 0.9(3) ^eV from the 
saturation of single-particle processes as detailed in Sec. IV B. With this the dimensionless Nambu spinor component 
V?/vq/S. that describes the spectral weight of the Shiba state at the impurity site remains the only unknown parameter 
in our model. Our results also predict the relation a+/3+/a_/3_ = 1 in the linear regime, which we can use as an 
additional check of the robustness of our theoretical description. According to the data shown in Fig. Sll(b) this 
relation is satisfied remarkably well over more than two decades of normal state conductance, throughout the linear 
regime. 

The theoretical conductance peak heights obtained numerically from Eqs. (S51) and (S52) are plotted in Fig. Sll(a) 
as a function of the normal state conductance 


Gn 


1 + 


(S89) 


While several features of the theoretical curves qualitatively agree with the experimental data in Fig. 3 of the main text 
there are also notable deviations. Most prominently, theory predicts a peak in a_ at intermediate tunneling strength 
absent in the experiment. We attribute this deviation to extrinsic broadening introduced by the measurement setup, 
e.g., due to radio frequency noise. Indeed, in the low coupling limit we find a peak width ic ~ 70 peV [see black curve in 
SlO(b)] exceeding Fi by more than an order of magnitude. To account for the broadening we convolute the theoretical 
dl/dV curves with a Gaussian of width w and plot the resulting peak heights in Fig. 3 of the main text. We find 
remarkable agreement with the experimental data over the entire range of normal state conductance and determine 
M^/r'oA ~ 0.23 from the fit. We associate deviations for /3_ at large normal state conductances ^ 0.02(2e^//i) 
with a multiple Andreev reflection resonance involving the second Shiba state at ei ~ 0.77 meV. A resonance occurs at 
eV = —(A-|-ei)/2 ~ —l.OOmeV and therefore overlaps with the thermal Shiba peak at eV = —(A —eo) — —1.13meV 
[see blue curve in Fig. SlO(b)]. 

In Fig. SI2(a), we plot the current measured at the position of the conductance peaks together with the theoretical 
curves. These fits use the same parameters as for the conductance hts, including the extrinsic broadening. We again 
hnd excellent agreement which corroborates that our model calculation correctly captures the essential tunneling 
processes. 
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normal state dl/dV (2e%) normal state d//d\/(2e%) 


Figure S12. (a) Measured current at 1.2 K at the two main Shiba peaks (eV^ = ±(A + eo)), and the corresponding two thermal 
peaks {eV = ±(A — eo))- The lines show the theoretical curves using the parameters mentioned in the text. The parameters 
are the same as for the ht of the differential conductance in Fig. 3 of the main text, (b) Measured current at 4.8 K. 


B. Quasiparticle lifetime and relaxation mechanism 

The quasiparticle (quasihole) lifetimes of the Shiba state are related to the inverse relaxation rates ti /2 = ^1^1/2- 
These rates could in principle be determined from the linewidth of the Shiba resonance at weak coupling (see Sec. Ill F). 
The measured linewidth, however, may be increased by an additional broadening from the measurement setup that 
sets the experimental energy resolution. The intrinsic linewidth of the data taken T = 1.2K is most likely well below 
the resolution of the experiment (see also Sec. IV A). 

A more robust way to determine the lifetime focuses on the strong-tunneling regime where the data remains 
unaffected by the energy resolution. As discussed in Sec. HIE, the single-particle current saturates when Ti < We 
assuming a value of = eri/ 2 / 3 /i at the main (thermal) thresholds. The current measured at the thermal threshold 
eV = A — eo shown in Fig. S12(b) indeed exhibits a plateau at strong tunneling. At eV = A — eo we extract 
a saturation current of 0.07(2) nA which yields relaxation rates r 2 = 0.9(3) peV and Fi ~ 4 F 2 = 4(1) peV. The 
corresponding lifetimes ti ~ 0.2ns and T 2 — 0.7ns are quoted in the main text. The current /_(A-eo) other 

thermal threshold does not saturate because of additional subgap features at strong tunneling discussed in Sec. IV A. 

Alternatively, we can determine Fi from the current at the main threshold. Because of the additional Andreev 
contribution, the current /_(A-i-eo) exhibits a shoulder instead of a plateau when the single-particle current saturates 
[cf. Fig. S9(b)]. The current at the shoulder is 0.4(2) nA and thus Fi = 5(3) peV, in agreement with the above value. 
While determining F 1 from the shoulder at the main thresholds typically has a larger uncertainty, it is furthermore 
subject to a systematic error due to extrinsic broadening introduced in the measurement as detailed in Sec. IV A. 
Note that the shoulder is absent in lA+to because the Andreev current dominates already when the single-particle 
current saturates. This behavior is well captured by a quantitative comparison between theory and experiment in 
Fig. S12(a). 

In order to gain insight into the relevant relaxation mechanisms in our system, we also study transport at a higher 
temperature of T = 4.8 K. While the qualitative temperature dependence of the relaxation rates ri /2 indicates 
thermal relaxation as discussed in the main text, a more quantitative analysis is required to assess the relevance 
of particular relaxation processes. At the higher temperature relaxation is strong enough that the single-particle 
current does not saturate in the tunneling regime Gm ^ 2e^//i as shown by the data in Fig. S12(b). We can instead 
determine Fi -|- r 2 from the peak width in the weak-coupling regime which we find to be ~ 0.16 meV. Assuming the 
extrinsic broadening is the same for both temperatures, we estimate Fi -|- r 2 ~ 0.2 meV [note that the unbroadened 
peak width is 0.7(Fi -|- F 2 )]. From the conductance peaks shown in the inset of Fig. 3 of the main text, we find 
a_//3+ = F 1 /F 2 ~ 1.6 in the linear regime and thus Fi ~ 120 peV which correspond to ti ~ 6ps. Hence we find the 
ratio of relaxation rates at the two temperatures Fi(4.8K)/Fi(1.2K) ~ 35. 

In Sec. IIF we have discussed quasiparticle transitions between Shiba state and quasiparticle continuum as possible 
relaxation mechanisms. Based on the corresponding relaxation rate in Eq. (S54), however, we would expect a ratio 
Fi(4.8K)/Fi(1.2K) ~ 10"^ (using eg = 0.20meV and A = 1.21meV at 4.8K), reflecting the exponential suppression 
of thermal relaxation at low temperatures. In view of the large discrepancy to the experimental value we can exclude 
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direct transitions to the quasiparticle continuum as the dominant relaxation process. An alternative scenario involves 
a cascade of transitions first to the next Shiba state at ei = 0.77 meV (ei = 0.68 meV at 4.8 K) and subsequent 
relaxation to the third Shiba state and the continuum. A rough estimate using Eq. (S54) with ei instead of A yields a 
ratio ri(4.8K)/ri(1.2K) ^ 10^. This result is much closer to the experimental value though a more detailed analysis 
of the population dynamics of the various Shiba states is required for a quantitative comparison. Note that a slightly 
higher temperature than 1.2 K could account for a significant part of this deviation. 


